/* Program to calculate PI */

#include <stdlib.h>             /* Standard Lib    ISOC  */
#include <stdarg.h>             /* Variable args   ISOC  */
#include <stdio.h>              /* I/O lib         ISOC  */
#include <gmp.h>                /* GNU GMP Library       */
#include <math.h>
#include <pthread.h>


typedef struct{
mpf_t  a_mp, b_mp, t_mp, x_mp, y_mp, temp_mp1,temp_mp2,temp_mp3,temp_mp4,temp_mp5, pi_mp;
}GaussLegendre;

void * passo1(void *arg){
	GaussLegendre *alg = arg;

	//a = (a + b) / 2;
	mpf_add (alg->a_mp, alg->a_mp, alg->b_mp);
	mpf_div_ui (alg->a_mp, alg->a_mp, 2);

	return NULL;
}



void * passo2(void *arg){
	GaussLegendre *alg = arg;


	//b = sqrt(b * y);
	mpf_mul (alg->temp_mp2 , alg->b_mp , alg->y_mp);
	mpf_sqrt (alg->b_mp, alg->temp_mp2);
	return NULL;
}


void * passo3(void *arg){
	GaussLegendre *alg = arg;


	//t = t - x*((y - a)*(y-a));
	mpf_sub (alg->temp_mp3, alg->y_mp, alg->a_mp);
	mpf_mul (alg->temp_mp3, alg->temp_mp3, alg->temp_mp3);
	mpf_mul (alg->temp_mp3, alg->x_mp, alg->temp_mp3);
	mpf_sub (alg->t_mp, alg->t_mp, alg->temp_mp3);
	return NULL;
}


void * passo4(void *arg){
	GaussLegendre *alg = arg;

	//x = 2*x;
	mpf_mul_ui (alg->x_mp, alg->x_mp, 2);
	return NULL;
}

void * passo5(void *arg){
	GaussLegendre *alg = arg;

	//y = a;
	mpf_set (alg->y_mp, alg->a_mp);
	return NULL;
}

int main()
{
//double a=1, b, t=0.25, x=1, y;
GaussLegendre alg;

pthread_t threads[5];

unsigned int its, i;





//b = 1/sqrt(2);

//printf("\n Entre com o numero de casas decimais de precisao: >");
//scanf("%ud", &its);

//seta a precisao do ponto flutuante
mpf_set_default_prec(1024 * 40000);
mpf_init_set_d(alg.a_mp,  1.0);

mpf_init_set_d(alg.t_mp,  0.25);
mpf_init_set_d(alg.x_mp,  0.5);
mpf_init_set_d(alg.y_mp,  0.0);
mpf_init_set_d(alg.temp_mp1,  0.0);
mpf_init_set_d(alg.temp_mp2,  0.0);
mpf_init_set_d(alg.temp_mp3,  0.0);
mpf_init_set_d(alg.temp_mp4,  0.0);
mpf_init_set_d(alg.temp_mp5,  0.0);
mpf_init_set_d(alg.pi_mp,  0.0);


mpf_sqrt_ui (alg.temp_mp1, 2);
mpf_ui_div (alg.temp_mp1, 1, alg.temp_mp1);
mpf_init_set(alg.b_mp, alg.temp_mp1);



for(i=0; i<22; i++)
{
	pthread_create(threads + 4, NULL, passo5, &alg);
	//y = a;
	
	pthread_create(threads + 3, NULL, passo4, &alg);
	//x = 2*x
	
	

	pthread_join( threads[4],NULL);
	pthread_create(threads , NULL, passo1, &alg);
	//a = (a + b) / 2;
	pthread_join( threads[0],NULL);
	pthread_create(threads + 1, NULL, passo2, &alg);
	//b = sqrt(b * y);

	pthread_join( threads[3],NULL);

	pthread_create(threads + 2, NULL, passo3, &alg);
	//t = t - x*((y - a)*(y-a));
	
	pthread_join( threads[1],NULL);
	
	pthread_join( threads[2],NULL);
	
	
}

//printf("\nPI = %.10000f\n\n", ((a+b)*(a+b))/(4*t));
mpf_add (alg.temp_mp1, alg.b_mp, alg.a_mp);
mpf_mul (alg.temp_mp1, alg.temp_mp1, alg.temp_mp1);
mpf_mul_ui (alg.pi_mp, alg.t_mp, 4);
mpf_div (alg.pi_mp, alg.temp_mp1, alg.pi_mp);


printf("\n");
mpf_out_str(NULL, 10, 10000000, alg.pi_mp);
printf("\n");
return(0);
}
